#Responsability:
#Calculate the likelihood

library(stats)

#TODO: check the for with more than one measurements in method2.R
Likelihood <- setRefClass(    
  "likelihood"
  
  , fields = list(
    output="data.frame",
    measurements="data.frame"
  )
  , methods = list(
    #
    #
    # Constructor
    #
    #
    initialize = function(..., 
                          output = "",
                          measurements = "")
    {
      #TODO: validate data 
      print("Initialize Likelihood")
      callSuper(...,
                output = output,
                measurements = measurements)
    },
    dummyRandom = function() {
      
      rnd <- runif(1, 1.0, 17.5)
      
      loginfo("Dummy Random")
      loginfo(rnd)
      
      rnd
    },
    #
    #
    # Calculate the logarithmic likelihood between output and measurements
    #
    #
    logLikelihood = function() {
      
      logdebug("METHOD IN: likelihood")
      
      
      # last day(s) 365 (and 366) ommitted because 7 days * 52 weeks are only 364 days
      # initialize weeklymean and -output to NaN to decide later, whether they are set or not 
      weeklymean   <- NaN
      weeklyoutput <- NaN
      weeklysd <- NaN
      
      #51 semanas? Se podria paralelizar
      for (i in 0:51) {
        logdebug(paste("likelihood iteration ", i))
        #show(colnames(measurements))
        # boolean vector with TRUE's in right days
        week_i <-  measurements$day >= 1 + i*7 & measurements$day <= 7 + i*7
        
        # if there are no measurements in one week the week is omitted
        if ( sum(week_i) > 0  ) {
          
          #Nacho Saca la media y el error de los datos de la semana dada (las semanas es un vector de booleanos)
          mean_standarderror <- function(measurementSubDF, week, type)
          {
            logdebug("SUBMETHOD IN: mean_standarderror")
            # mean of week is weighted by number of daily replicates V1=year, V2=day, V3=measurement, V4=sd, V5=no of replicates
            
            totalreplicates <- sum( measurementSubDF$noRep[ week ] )
            weekmean <- sum( measurementSubDF$noRep[ week ]  * measurementSubDF$n_n2o[ week ] )  / totalreplicates
            
            # sigma := sd / sqrt n, if only 1 replicate the value is used as sd
            if ( totalreplicates > 1 ){
              betweenDAYvariance <- sum( ( measurementSubDF$n_n2o[ week ] - weekmean )^2 * measurementSubDF$noRep[ week ] )
              
              
              withinDayvariance  <- sum( ( measurementSubDF$noRep[ week ]-1) * (measurementSubDF$n_n2o_std[week])^2   )
              totalvariance <- withinDayvariance + betweenDAYvariance
            } else {
              totalvariance <- measurementSubDF$n_n2o_std
              
            }
            weekerror <- sqrt( totalvariance / totalreplicates )
            
            
            logdebug("SUBMETHOD OUT: mean_standarderror")
            
            c(weekmean, weekerror)
          }
          
          # week_i es un vector de booleanos para seleccionar los datos de la semana
          #measurementSubDF son todos los datos
          results <- mean_standarderror(measurements, week_i)
          #show(results)
          #Nacho se van insertando los datos calculados en el vector weeklymean y weeklysd
          if (is.nan( weeklymean[1] ))
          {
            weeklymean   <- results[1]
            weeklysd     <- results[2]
          }
          else { 
            weeklymean   <- c(weeklymean, results[1])
            weeklysd     <- c(weeklysd, results[2]) 
          }
          
          if (is.nan( weeklyoutput[1] )) {
            weeklyoutput <- mean(output[week_i])
          } else {
            weeklyoutput <- c(weeklyoutput, mean(output[week_i]))
          }
          
        } # endif sum(week_i)
        
      }
      
      #¿Que tenemos hasta ahora?
      # - Weeklymean: es un vector que contiene las medias de cada semana de los measurements
      # - Weeklysd: es un vector que contiene los sd de cada semana de los measurements
      # - weeklyoutput: los mismo, un vector que contiene la media del output para cada semana.
      
      #Nacho el nombre loglikelihood2 entiendo que es el logaritmo de la probabilidad al cuadrado, pero no estoy seguro
      
      likelihood <- sum(loglikelihood2( weeklyoutput, weeklymean , weeklysd))
      
      #a es la suma de todas las probabilidades p(D|Teta) resultantes de cada semana
      
      
      if ( is.infinite(likelihood) )  likelihood <- -9999999
      if ( is.na(likelihood)) likelihood <- -8888888
      if ( is.nan(likelihood)) likelihood <- -7777777
      
      logdebug("METHOD OUT: likelihood")
      #      show(likelihood)
      
      likelihood
      
    },  
    # loglikelihood with normal distribution (not sivia function)
    #- modelOutput: outputs of the models with weekly means
    #- meas: measurements: column with data in the measurement file (ex: co2)
    #- measurementsSD: measurements Standard Deviation: column with the sd in the measurement file.
    loglikelihood2 = function(modelOutput, meas, measurementsSD) {
      
      loginfo("METHOD IN: loglikelihood2")
      
      #        show(measurementsSD)
      
      #Nacho Mod sum(meaSD <= 0) now is sum(meaSD < 0)
      if ((length(measurementsSD) < 1 ) | ( sum(measurementsSD < 0) > 0 )) {
        stop("Standard Deviation of measurement (meaSD) is not positive (or does not exist)"  )
      }
      
      if (length(meas) == 1 & length(measurementsSD) == 1 & length(modelOutput) != 1) {
        meas   <- vector(length= length(modelOutput), mode="numeric") + meas
        measurementsSD <- vector(length= length(modelOutput), mode="numeric") + measurementsSD
      }
      
      if (length(meas) != length(measurementsSD)) {
        stop("Error: Length of mea and meaSD are not equal")
      } else {
        
        likeli <- dnorm(modelOutput, meas, measurementsSD) 
        
        logLi <- log(likeli)
        
        loginfo("METHOD OUT: loglikelihood2")
        
        logLi
      }
      
    }
    
    
    
    
    
  )#End methods List 
  
)#End RefClass